rm(list = ls())
library(tidyverse)
library(kableExtra)
library(ggplot2)


load("data/results_pooled_est.RData")

load("data/results_var_comp.RData")

1 Table

table <- perf_crit  |>
  group_by(model, variance, var_combo_graph) |>
  summarize(
    mean_pb = mean(bias_icc),
    var_pb = var(bias_icc),
    min_pb = min(bias_icc),
    max_pb = max(bias_icc),
    
    # mean_pb_mcse = mean(bias_icc_mcse),
    # var_pb_mcse = var(bias_icc_mcse),
    # min_pb_mcse = min(bias_icc_mcse),
    # max_pb_mcse = max(bias_icc_mcse),
    
    
    mean_rpb = mean(rel_bias_icc),
    var_rpb = sd(rel_bias_icc),
    min_rpb = min(rel_bias_icc),
    max_rpb = max(rel_bias_icc),
    
    mean_rmse = mean(RMSE),
    var_rmse = sd(RMSE),
    min_rmse = min(RMSE),
    max_rmse = max(RMSE),
    
    mean_mse = mean(MSE),
    var_mse = sd(MSE),
    min_mse = min(MSE),
    max_mse = max(MSE),
    
    mean_rseb = mean(rel_bias_se_ICC),
    var_rseb = sd(rel_bias_se_ICC),
    min_rseb = min(rel_bias_se_ICC),
    max_rseb = max(rel_bias_se_ICC),
    
    # mean_rpb_mcse = mean(rel_bias_icc_mcse),
    # var_rpb_mcse = var(rel_bias_icc_mcse),
    # min_rpb_mcse = min(rel_bias_icc_mcse),
    # max_rpb_mcse = max(rel_bias_icc_mcse
    .groups = "drop"
  ) |>
  arrange(var_combo_graph) |>
  rename(
    Estimation = model,
    "Variance Formulae" = variance,
    "ICC Parameter" = var_combo_graph,
    "Parameter Bias - mean" = mean_pb,
    "Parameter Bias - variance" = var_pb,
    "Parameter Bias - min." = min_pb,
    "Parameter Bias - max" = max_pb,
    "Relative Parameter Bias - mean" = mean_rpb,
    "Relative Parameter Bias - variance" = var_rpb,
    "Relative Parameter Bias - min." = min_rpb,
    "Relative Parameter Bias - max" = max_rpb,
    "Root Mean Squared Error - mean" = mean_rmse,
    "Root Mean Squared Error - variance" = var_rmse,
    "Root Mean Squared Error - min." = min_rmse,
    "Root Mean Squared Error - max" = max_rmse,
    "Mean Squared Error - mean" = mean_mse,
    "Mean Squared Error - variance" = var_mse,
    "Mean Squared Error - min." = min_mse,
    "Mean Squared Error - max" = max_mse,
    "Relative Standard Error Bias - mean" = mean_rseb,
    "Relative Standard Error Bias - variance" = var_rseb,
    "Relative Standard Error Bias - min." = min_rseb,
    "Relative Standard Error Bias - max" = max_rseb#,
    
    #
    # mean_rpb_ratio = mean(rel_bias_icc_alt),
    # min_rpb_ratio = min(rel_bias_icc_alt),
    # max_rpb_ratio = max(rel_bias_icc_alt),
    # sd_rpb_ratio = sd(rel_bias_icc_alt)
    
  )

  write.csv(table, file = "perform_crit_summary.csv")
  
  table |> 
   kbl() |>
   kable_styling() |> 
   scroll_box(width = "900px", height = "500px")
Estimation Variance Formulae ICC Parameter Parameter Bias - mean Parameter Bias - variance Parameter Bias - min. Parameter Bias - max Relative Parameter Bias - mean Relative Parameter Bias - variance Relative Parameter Bias - min. Relative Parameter Bias - max Root Mean Squared Error - mean Root Mean Squared Error - variance Root Mean Squared Error - min. Root Mean Squared Error - max Mean Squared Error - mean Mean Squared Error - variance Mean Squared Error - min. Mean Squared Error - max Relative Standard Error Bias - mean Relative Standard Error Bias - variance Relative Standard Error Bias - min. Relative Standard Error Bias - max
REML Donner 0.05 -0.0055206 3.10e-06 -0.0086502 -0.0033179 -0.1104120 0.0350042 -0.1730048 -0.0663589 0.0068103 0.0018541 0.0038674 0.0109862 0.0000497 0.0000262 1.50e-05 0.0001207 -0.0418569 0.0298111 -0.0980127 0.0093520
REML Fisher 0.05 -0.0055134 3.10e-06 -0.0084706 -0.0032662 -0.1102687 0.0351356 -0.1694117 -0.0653237 0.0068032 0.0018560 0.0038551 0.0108288 0.0000497 0.0000261 1.49e-05 0.0001173 -0.0419004 0.0297949 -0.0973853 0.0089611
REML Fisher TF 0.05 -0.0038014 1.30e-06 -0.0055870 -0.0020606 -0.0760285 0.0224458 -0.1117401 -0.0412123 0.0053695 0.0015421 0.0026996 0.0083844 0.0000312 0.0000170 7.30e-06 0.0000703 -0.0021009 0.0235978 -0.0595003 0.0393137
REML Hedges 0.05 -0.0054133 3.00e-06 -0.0084060 -0.0032763 -0.1082661 0.0343600 -0.1681204 -0.0655256 0.0067228 0.0018241 0.0038339 0.0108028 0.0000485 0.0000253 1.47e-05 0.0001167 -0.0406881 0.0295053 -0.0949330 0.0108128
REML Smith 0.05 -0.0056378 3.20e-06 -0.0088097 -0.0033974 -0.1127560 0.0358532 -0.1761949 -0.0679483 0.0069050 0.0018834 0.0039056 0.0110766 0.0000512 0.0000269 1.53e-05 0.0001227 -0.0428401 0.0300361 -0.0980784 0.0085310
REML Swiger 0.05 -0.0055107 3.10e-06 -0.0084615 -0.0032660 -0.1102139 0.0350982 -0.1692307 -0.0653193 0.0068010 0.0018545 0.0038549 0.0108251 0.0000496 0.0000261 1.49e-05 0.0001172 -0.0418829 0.0297847 -0.0971925 0.0089368
RVE Donner 0.05 -0.0055777 3.00e-06 -0.0085649 -0.0033146 -0.1115539 0.0346341 -0.1712981 -0.0662929 0.0068763 0.0018542 0.0038117 0.0109456 0.0000507 0.0000262 1.45e-05 0.0001198 -0.0541107 0.0453614 -0.1367120 0.0244123
RVE Fisher 0.05 -0.0055726 3.00e-06 -0.0083736 -0.0033254 -0.1114524 0.0347894 -0.1674713 -0.0665089 0.0068696 0.0018546 0.0037992 0.0107733 0.0000506 0.0000261 1.44e-05 0.0001161 -0.0532168 0.0450808 -0.1369035 0.0266975
RVE Fisher TF 0.05 -0.0038014 1.30e-06 -0.0055868 -0.0020606 -0.0760289 0.0224448 -0.1117352 -0.0412113 0.0053696 0.0015422 0.0026996 0.0083850 0.0000312 0.0000170 7.30e-06 0.0000703 -0.0020099 0.0235566 -0.0594670 0.0397569
RVE Hedges 0.05 -0.0054668 2.90e-06 -0.0083203 -0.0032798 -0.1093365 0.0339513 -0.1664052 -0.0655957 0.0067853 0.0018239 0.0037795 0.0107628 0.0000493 0.0000254 1.43e-05 0.0001158 -0.0540002 0.0443273 -0.1361919 0.0245538
RVE Smith 0.05 -0.0056991 3.20e-06 -0.0087006 -0.0033549 -0.1139811 0.0355106 -0.1740112 -0.0670980 0.0069730 0.0018801 0.0038499 0.0110027 0.0000521 0.0000269 1.48e-05 0.0001211 -0.0530666 0.0462642 -0.1369892 0.0268108
RVE Swiger 0.05 -0.0055698 3.00e-06 -0.0083648 -0.0033253 -0.1113955 0.0347507 -0.1672966 -0.0665057 0.0068672 0.0018531 0.0037991 0.0107691 0.0000505 0.0000261 1.44e-05 0.0001160 -0.0532047 0.0450506 -0.1368946 0.0266430
REML Donner 0.10 -0.0070980 4.50e-06 -0.0103775 -0.0043971 -0.0709803 0.0211630 -0.1037751 -0.0439705 0.0086478 0.0022452 0.0050989 0.0130228 0.0000797 0.0000395 2.60e-05 0.0001696 -0.0277921 0.0230186 -0.0795846 0.0331137
REML Fisher 0.10 -0.0071406 4.60e-06 -0.0103323 -0.0044362 -0.0714056 0.0214930 -0.1033234 -0.0443619 0.0086834 0.0022679 0.0051412 0.0131500 0.0000804 0.0000401 2.64e-05 0.0001729 -0.0281511 0.0229956 -0.0802938 0.0325047
REML Fisher TF 0.10 -0.0049637 1.70e-06 -0.0073567 -0.0025374 -0.0496370 0.0130568 -0.0735668 -0.0253741 0.0068912 0.0017562 0.0034347 0.0109561 0.0000505 0.0000250 1.18e-05 0.0001200 -0.0026267 0.0205246 -0.0413956 0.0444442
REML Hedges 0.10 -0.0069785 4.40e-06 -0.0101173 -0.0043557 -0.0697847 0.0209693 -0.1011728 -0.0435567 0.0085501 0.0022251 0.0050570 0.0128375 0.0000780 0.0000386 2.56e-05 0.0001648 -0.0270140 0.0231271 -0.0794918 0.0341317
REML Smith 0.10 -0.0073212 4.80e-06 -0.0108172 -0.0044442 -0.0732116 0.0220040 -0.1081715 -0.0444417 0.0088344 0.0023074 0.0051495 0.0131570 0.0000833 0.0000415 2.65e-05 0.0001731 -0.0288332 0.0233585 -0.0802933 0.0325404
REML Swiger 0.10 -0.0071367 4.60e-06 -0.0103242 -0.0044352 -0.0713670 0.0214699 -0.1032416 -0.0443520 0.0086802 0.0022661 0.0051401 0.0131430 0.0000804 0.0000400 2.64e-05 0.0001727 -0.0281266 0.0229980 -0.0802703 0.0325331
RVE Donner 0.10 -0.0070986 4.50e-06 -0.0102447 -0.0043907 -0.0709858 0.0212475 -0.1024472 -0.0439069 0.0086500 0.0022356 0.0051176 0.0129352 0.0000797 0.0000393 2.62e-05 0.0001673 -0.0126710 0.0331201 -0.0803805 0.0558477
RVE Fisher 0.10 -0.0071480 4.70e-06 -0.0102092 -0.0044334 -0.0714798 0.0216354 -0.1020923 -0.0443341 0.0086913 0.0022620 0.0051666 0.0130718 0.0000805 0.0000398 2.67e-05 0.0001709 -0.0123514 0.0335987 -0.0793201 0.0568207
RVE Fisher TF 0.10 -0.0049635 1.70e-06 -0.0073566 -0.0025376 -0.0496354 0.0130562 -0.0735663 -0.0253758 0.0068911 0.0017563 0.0034347 0.0109561 0.0000505 0.0000250 1.18e-05 0.0001200 -0.0030384 0.0205266 -0.0416797 0.0444772
RVE Hedges 0.10 -0.0069702 4.40e-06 -0.0099708 -0.0043473 -0.0697021 0.0209913 -0.0997078 -0.0434734 0.0085450 0.0022116 0.0050710 0.0127432 0.0000778 0.0000382 2.57e-05 0.0001624 -0.0137747 0.0322466 -0.0801877 0.0542139
RVE Smith 0.10 -0.0073410 4.90e-06 -0.0107213 -0.0044420 -0.0734098 0.0222217 -0.1072128 -0.0444202 0.0088525 0.0023053 0.0051759 0.0130782 0.0000836 0.0000413 2.68e-05 0.0001710 -0.0106762 0.0344995 -0.0792582 0.0655423
RVE Swiger 0.10 -0.0071438 4.70e-06 -0.0102015 -0.0044323 -0.0714385 0.0216096 -0.1020149 -0.0443229 0.0086878 0.0022600 0.0051653 0.0130641 0.0000805 0.0000398 2.67e-05 0.0001707 -0.0123816 0.0335700 -0.0793435 0.0567584
REML Donner 0.15 -0.0081131 5.70e-06 -0.0114168 -0.0051556 -0.0540873 0.0159411 -0.0761119 -0.0343704 0.0100361 0.0025978 0.0061445 0.0152725 0.0001073 0.0000531 3.78e-05 0.0002332 -0.0244915 0.0239332 -0.0795617 0.0473734
REML Fisher 0.15 -0.0081880 5.90e-06 -0.0116082 -0.0051958 -0.0545870 0.0162508 -0.0773881 -0.0346387 0.0100986 0.0026305 0.0061904 0.0152935 0.0001088 0.0000540 3.83e-05 0.0002339 -0.0248300 0.0240200 -0.0803382 0.0469676
REML Fisher TF 0.15 -0.0063898 3.70e-06 -0.0095806 -0.0032772 -0.0425987 0.0127480 -0.0638704 -0.0218479 0.0086216 0.0022393 0.0043784 0.0137339 0.0000792 0.0000399 1.92e-05 0.0001886 -0.0077120 0.0244446 -0.0627136 0.0337777
REML Hedges 0.15 -0.0079890 5.70e-06 -0.0111627 -0.0051144 -0.0532603 0.0158501 -0.0744181 -0.0340960 0.0099354 0.0025843 0.0060984 0.0150839 0.0001053 0.0000521 3.72e-05 0.0002275 -0.0237531 0.0239181 -0.0797239 0.0472985
REML Smith 0.15 -0.0083992 6.20e-06 -0.0119814 -0.0052049 -0.0559948 0.0166630 -0.0798761 -0.0346993 0.0102747 0.0026759 0.0062003 0.0156848 0.0001126 0.0000559 3.84e-05 0.0002460 -0.0255192 0.0242548 -0.0841158 0.0469529
REML Swiger 0.15 -0.0081836 5.90e-06 -0.0115976 -0.0051948 -0.0545571 0.0162335 -0.0773176 -0.0346319 0.0100949 0.0026284 0.0061892 0.0152860 0.0001087 0.0000539 3.83e-05 0.0002337 -0.0248006 0.0240173 -0.0803106 0.0469787
RVE Donner 0.15 -0.0081877 6.00e-06 -0.0115896 -0.0051975 -0.0545850 0.0163212 -0.0772643 -0.0346503 0.0100978 0.0026211 0.0062413 0.0152677 0.0001087 0.0000537 3.90e-05 0.0002331 -0.0006777 0.0294194 -0.0746190 0.0730476
RVE Fisher 0.15 -0.0082726 6.30e-06 -0.0116808 -0.0052429 -0.0551507 0.0166864 -0.0778721 -0.0349526 0.0101690 0.0026597 0.0062968 0.0152887 0.0001103 0.0000547 3.97e-05 0.0002337 -0.0001122 0.0298060 -0.0747036 0.0738986
RVE Fisher TF 0.15 -0.0063894 3.70e-06 -0.0095800 -0.0032776 -0.0425958 0.0127463 -0.0638666 -0.0218505 0.0086213 0.0022391 0.0043784 0.0137329 0.0000792 0.0000399 1.92e-05 0.0001886 -0.0082602 0.0244407 -0.0634913 0.0330649
RVE Hedges 0.15 -0.0080501 5.90e-06 -0.0112830 -0.0051529 -0.0536675 0.0161842 -0.0752201 -0.0343530 0.0099856 0.0026029 0.0061880 0.0150658 0.0001063 0.0000526 3.83e-05 0.0002270 -0.0017014 0.0290064 -0.0748050 0.0720901
RVE Smith 0.15 -0.0085070 6.70e-06 -0.0122454 -0.0052528 -0.0567136 0.0171966 -0.0816357 -0.0350185 0.0103658 0.0027158 0.0063082 0.0157074 0.0001147 0.0000570 3.98e-05 0.0002467 0.0017423 0.0305753 -0.0730971 0.0740775
RVE Swiger 0.15 -0.0082676 6.20e-06 -0.0116683 -0.0052417 -0.0551172 0.0166663 -0.0777885 -0.0349448 0.0101648 0.0026574 0.0062954 0.0152810 0.0001102 0.0000547 3.96e-05 0.0002335 -0.0001479 0.0297885 -0.0747047 0.0738751
REML Donner 0.25 -0.0086569 6.50e-06 -0.0117173 -0.0055425 -0.0346275 0.0102083 -0.0468692 -0.0221701 0.0114385 0.0029811 0.0068325 0.0168380 0.0001395 0.0000693 4.67e-05 0.0002835 -0.0133702 0.0258710 -0.0560011 0.0644447
REML Fisher 0.25 -0.0087611 6.80e-06 -0.0118725 -0.0055795 -0.0350445 0.0104341 -0.0474900 -0.0223181 0.0115223 0.0030168 0.0068786 0.0169508 0.0001417 0.0000705 4.73e-05 0.0002873 -0.0137947 0.0258242 -0.0560912 0.0642463
REML Fisher TF 0.25 -0.0089397 9.20e-06 -0.0134210 -0.0044222 -0.0357587 0.0121512 -0.0536839 -0.0176887 0.0116512 0.0031704 0.0057559 0.0177719 0.0001456 0.0000743 3.31e-05 0.0003158 0.0010399 0.0253765 -0.0712062 0.0528747
REML Hedges 0.25 -0.0085434 6.50e-06 -0.0114737 -0.0054110 -0.0341736 0.0101869 -0.0458949 -0.0216439 0.0113503 0.0029770 0.0067886 0.0166920 0.0001375 0.0000686 4.61e-05 0.0002786 -0.0127728 0.0258112 -0.0561434 0.0640485
REML Smith 0.25 -0.0089778 7.10e-06 -0.0124868 -0.0055883 -0.0359112 0.0106921 -0.0499473 -0.0223533 0.0116977 0.0030550 0.0071475 0.0172391 0.0001460 0.0000724 5.11e-05 0.0002972 -0.0146388 0.0259807 -0.0562466 0.0648744
REML Swiger 0.25 -0.0087564 6.80e-06 -0.0118605 -0.0055775 -0.0350257 0.0104238 -0.0474420 -0.0223099 0.0115185 0.0030150 0.0068773 0.0169447 0.0001416 0.0000705 4.73e-05 0.0002871 -0.0137662 0.0258247 -0.0560896 0.0642467
RVE Donner 0.25 -0.0088005 7.00e-06 -0.0119670 -0.0055898 -0.0352018 0.0105707 -0.0478682 -0.0223592 0.0115555 0.0030182 0.0069429 0.0168908 0.0001425 0.0000705 4.82e-05 0.0002853 0.0127902 0.0302545 -0.0419399 0.0911629
RVE Fisher 0.25 -0.0089165 7.30e-06 -0.0121348 -0.0056306 -0.0356660 0.0108344 -0.0485391 -0.0225224 0.0116499 0.0030605 0.0069931 0.0170182 0.0001449 0.0000720 4.89e-05 0.0002896 0.0132855 0.0303535 -0.0416434 0.0917111
RVE Fisher TF 0.25 -0.0089388 9.20e-06 -0.0134198 -0.0044227 -0.0357552 0.0121479 -0.0536793 -0.0176909 0.0116504 0.0031698 0.0057559 0.0177705 0.0001456 0.0000743 3.31e-05 0.0003158 0.0003851 0.0253995 -0.0723426 0.0523955
RVE Hedges 0.25 -0.0086741 6.90e-06 -0.0117676 -0.0054521 -0.0346965 0.0105305 -0.0470704 -0.0218084 0.0114562 0.0030115 0.0068941 0.0167275 0.0001401 0.0000698 4.75e-05 0.0002798 0.0122270 0.0302147 -0.0425272 0.0901897
RVE Smith 0.25 -0.0091580 7.80e-06 -0.0128584 -0.0056401 -0.0366319 0.0111526 -0.0514337 -0.0225603 0.0118476 0.0031086 0.0072980 0.0173469 0.0001498 0.0000743 5.33e-05 0.0003009 0.0143791 0.0307191 -0.0382305 0.0957725
RVE Swiger 0.25 -0.0089112 7.30e-06 -0.0121212 -0.0056289 -0.0356449 0.0108223 -0.0484848 -0.0225155 0.0116456 0.0030585 0.0069917 0.0170112 0.0001448 0.0000719 4.89e-05 0.0002894 0.0132648 0.0303513 -0.0416562 0.0916901
REML Donner 0.50 -0.0046887 2.60e-06 -0.0070313 -0.0028291 -0.0093774 0.0032041 -0.0140626 -0.0056582 0.0101298 0.0032817 0.0055367 0.0166666 0.0001132 0.0000722 3.07e-05 0.0002778 -0.0225623 0.0203117 -0.0637735 0.0272628
REML Fisher 0.50 -0.0047143 2.60e-06 -0.0070271 -0.0028216 -0.0094286 0.0032176 -0.0140541 -0.0056432 0.0101465 0.0032853 0.0055342 0.0166760 0.0001135 0.0000724 3.06e-05 0.0002781 -0.0230695 0.0203648 -0.0647129 0.0267452
REML Fisher TF 0.50 -0.0119975 2.09e-05 -0.0183252 -0.0056447 -0.0239950 0.0091452 -0.0366503 -0.0112893 0.0152137 0.0044832 0.0076579 0.0232952 0.0002511 0.0001361 5.86e-05 0.0005427 -0.0005446 0.0220473 -0.0385554 0.0643555
REML Hedges 0.50 -0.0046951 2.60e-06 -0.0070255 -0.0028390 -0.0093901 0.0031966 -0.0140509 -0.0056780 0.0101283 0.0032786 0.0055423 0.0166487 0.0001131 0.0000721 3.07e-05 0.0002772 -0.0221176 0.0202850 -0.0622546 0.0274322
REML Smith 0.50 -0.0047318 2.60e-06 -0.0070290 -0.0028227 -0.0094637 0.0032358 -0.0140579 -0.0056454 0.0101635 0.0032911 0.0055352 0.0166775 0.0001139 0.0000726 3.06e-05 0.0002781 -0.0237785 0.0204854 -0.0647166 0.0266716
REML Swiger 0.50 -0.0047139 2.60e-06 -0.0070271 -0.0028220 -0.0094278 0.0032170 -0.0140542 -0.0056440 0.0101461 0.0032851 0.0055343 0.0166753 0.0001135 0.0000724 3.06e-05 0.0002781 -0.0230428 0.0203606 -0.0646352 0.0267591
RVE Donner 0.50 -0.0046818 2.60e-06 -0.0070300 -0.0028209 -0.0093635 0.0031993 -0.0140600 -0.0056417 0.0101308 0.0032846 0.0055340 0.0166750 0.0001132 0.0000723 3.06e-05 0.0002781 -0.0074028 0.0190845 -0.0471985 0.0403454
RVE Fisher 0.50 -0.0047077 2.60e-06 -0.0070259 -0.0028127 -0.0094154 0.0032134 -0.0140518 -0.0056254 0.0101480 0.0032885 0.0055313 0.0166846 0.0001136 0.0000725 3.06e-05 0.0002784 -0.0073879 0.0190841 -0.0475849 0.0404717
RVE Fisher TF 0.50 -0.0119960 2.09e-05 -0.0183229 -0.0056450 -0.0239920 0.0091421 -0.0366457 -0.0112900 0.0152124 0.0044821 0.0076581 0.0232928 0.0002511 0.0001360 5.86e-05 0.0005426 -0.0012317 0.0220094 -0.0396525 0.0637564
RVE Hedges 0.50 -0.0046883 2.50e-06 -0.0070241 -0.0028373 -0.0093766 0.0031917 -0.0140481 -0.0056747 0.0101293 0.0032814 0.0055399 0.0166566 0.0001131 0.0000722 3.07e-05 0.0002774 -0.0073436 0.0190830 -0.0461815 0.0404359
RVE Smith 0.50 -0.0047256 2.60e-06 -0.0070279 -0.0028138 -0.0094511 0.0032325 -0.0140558 -0.0056276 0.0101656 0.0032946 0.0055323 0.0166861 0.0001140 0.0000727 3.06e-05 0.0002784 -0.0073035 0.0191023 -0.0475624 0.0404760
RVE Swiger 0.50 -0.0047073 2.60e-06 -0.0070260 -0.0028132 -0.0094145 0.0032127 -0.0140519 -0.0056264 0.0101476 0.0032883 0.0055314 0.0166839 0.0001136 0.0000725 3.06e-05 0.0002784 -0.0073856 0.0190835 -0.0475400 0.0404697
REML Donner 0.90 0.0026873 8.00e-07 0.0015520 0.0042820 0.0029859 0.0009654 0.0017245 0.0047578 0.0049811 0.0014917 0.0026439 0.0083876 0.0000270 0.0000161 7.00e-06 0.0000704 -0.0482666 0.0255422 -0.1002644 0.0095753
REML Fisher 0.90 0.0027883 9.00e-07 0.0015924 0.0044530 0.0030982 0.0010255 0.0017693 0.0049478 0.0050419 0.0015024 0.0026784 0.0084627 0.0000276 0.0000163 7.20e-06 0.0000716 -0.0489632 0.0256952 -0.1015855 0.0088566
REML Fisher TF 0.90 -0.0029852 4.00e-06 -0.0064162 0.0003189 -0.0033169 0.0022094 -0.0071291 0.0003544 0.0054809 0.0017145 0.0024727 0.0088977 0.0000329 0.0000193 6.10e-06 0.0000792 -0.0077060 0.0240620 -0.0603324 0.0479065
REML Hedges 0.90 0.0025850 7.00e-07 0.0015133 0.0040539 0.0028723 0.0009534 0.0016814 0.0045043 0.0049245 0.0014962 0.0026133 0.0082892 0.0000264 0.0000160 6.80e-06 0.0000687 -0.0477363 0.0256486 -0.0983524 0.0109008
REML Smith 0.90 0.0029853 1.00e-06 0.0017077 0.0050108 0.0033170 0.0011057 0.0018974 0.0055676 0.0051620 0.0015152 0.0026858 0.0087250 0.0000289 0.0000168 7.20e-06 0.0000761 -0.0501228 0.0258765 -0.1054508 0.0077634
REML Swiger 0.90 0.0027842 8.00e-07 0.0015912 0.0044430 0.0030935 0.0010227 0.0017680 0.0049367 0.0050394 0.0015019 0.0026775 0.0084583 0.0000276 0.0000163 7.20e-06 0.0000715 -0.0489314 0.0256892 -0.1015095 0.0088963
RVE Donner 0.90 0.0026360 7.00e-07 0.0015095 0.0041099 0.0029289 0.0009001 0.0016772 0.0045665 0.0049578 0.0014872 0.0026054 0.0083483 0.0000267 0.0000160 6.80e-06 0.0000697 -0.0621029 0.0269319 -0.1215069 0.0034897
RVE Fisher 0.90 0.0027390 7.00e-07 0.0015517 0.0042677 0.0030433 0.0009601 0.0017241 0.0047419 0.0050183 0.0014925 0.0026413 0.0084130 0.0000274 0.0000162 7.00e-06 0.0000708 -0.0612664 0.0273264 -0.1211787 0.0040406
RVE Fisher TF 0.90 -0.0029847 4.00e-06 -0.0064147 0.0003188 -0.0033164 0.0022088 -0.0071275 0.0003542 0.0054806 0.0017143 0.0024727 0.0088978 0.0000329 0.0000193 6.10e-06 0.0000792 -0.0082505 0.0240534 -0.0601349 0.0468299
RVE Hedges 0.90 0.0025347 6.00e-07 0.0014704 0.0039010 0.0028163 0.0008872 0.0016338 0.0043345 0.0049028 0.0014949 0.0025746 0.0082672 0.0000262 0.0000159 6.60e-06 0.0000683 -0.0626932 0.0265652 -0.1216705 0.0031726
RVE Smith 0.90 0.0029397 9.00e-07 0.0017278 0.0047894 0.0032663 0.0010429 0.0019198 0.0053215 0.0051383 0.0014967 0.0026489 0.0086381 0.0000286 0.0000165 7.00e-06 0.0000746 -0.0596093 0.0284146 -0.1197092 0.0041511
RVE Swiger 0.90 0.0027347 7.00e-07 0.0015506 0.0042585 0.0030386 0.0009573 0.0017229 0.0047317 0.0050157 0.0014922 0.0026403 0.0084092 0.0000273 0.0000162 7.00e-06 0.0000707 -0.0613019 0.0273074 -0.1212003 0.0040281

2 Pooled Estimate

2.1 Interactions

2.1.1 True ICC vs num pooled ICCs

perf_crit  |>  
  ggplot(aes(x = variance, y = rel_bias_icc_median, fill = model)) +
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_numICC.png", width = 6.8, height = 5, units = "in")
perf_crit  |>  
  ggplot(aes(x = variance, y = rel_bias_icc_alt, fill = model)) + 
   geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RPB - ratio", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_ratio_numICC.png", width = 6.8, height = 5, units = "in")
# perf_crit_MCSE <- perf_crit |> 
#   group_by(var_combo_graph, icc_est_n) |>  
#   summarise(mcse = max(bias_icc_mcse))


perf_crit  |>  
  ggplot(aes(x = variance, y = bias_icc_median, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid", colour = "black")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_numICC.png", width = 6.8, height = 5, units = "in")


# ggplot(perf_crit, aes(var_combo_graph_f, bias_icc, col=model ) ) +
#   facet_grid(variance ~ icc_est_n, 
#              labeller = label_bquote(rows = .(variance), 
#                                             cols = k ==.(icc_est_n))) +
#   geom_point( alpha=0.10,
#               position = position_dodge(width = 0.04) ) +
#   geom_smooth( se=FALSE ) + 
#   geom_hline( yintercept = 0 ) +
#   theme_minimal()


# ggplot(perf_crit, aes(var_combo_graph_f, bias_icc, col=skewness ) ) +
#   facet_grid(variance ~ icc_est_n, 
#              labeller = label_bquote(rows = .(variance), 
#                                             cols = k ==.(icc_est_n))) +
#   geom_point( alpha=0.10,
#               position = position_dodge(width = 0.04) ) +
#   geom_smooth( se=FALSE ) + 
#   geom_hline( yintercept = 0 ) +
#   theme_minimal()


perf_crit  |>  
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_numICC.png", width = 6.8, height = 5, units = "in")

2.1.2 True ICC vs tau (between study heterogeneity)

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "dashed")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_tau.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_tau.png", width = 6.8, height = 5, units = "in")



perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed") + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_tau.png", width = 6.8, height = 5, units = "in")

2.1.3 True ICC vs number of clusters

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "dashed")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = j== .(nj_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_j.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid")  + 
  # geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = j== .(nj_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_j.png", width = 6.8, height = 5, units = "in")


perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = j== .(nj_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_j.png", width = 6.8, height = 5, units = "in")

2.1.4 True ICC vs cluster size

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "solid")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
   geom_hline(yintercept = 0, linetype = "solid")  + 
  # geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_nj.png", width = 6.8, height = 5, units = "in")


perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RMSE_nj.png", width = 6.8, height = 5, units = "in")

2.1.5 True ICC vs Degree of unbalance

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
   geom_hline(yintercept = .05, linetype = "dashed")  + 
   geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RPB_n_bar_prop.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  # geom_hline(yintercept = .05, linetype = "dashed")  + 
  # geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("PB_n_bar_prop.png", width = 6.8, height = 5, units = "in")


perf_crit  |> 
  ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols =zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RMSE_n_bar_prop.png", width = 6.8, height = 5, units = "in")

2.1.6 number clusters vs cluster size

This is how the variance formulae perform across all generating ICC values.

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RPB_j_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("PB_j_nj.png", width = 6.8, height = 5, units = "in")

perf_crit  |> ggplot(aes(x = variance, y = RMSE, fill = model)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RMSE_j_nj.png", width = 6.8, height = 5, units = "in")
perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_j_nj_nomodel.png", width = 6.8, height = 5, units = "in")

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("PB_j_nj_nomodel.png", width = 6.8, height = 5, units = "in")

perf_crit  |> ggplot(aes(x = variance, y = RMSE, fill = nj_size)) + 
  geom_hline(yintercept = .00, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RMSE", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RMSE_j_nj.png", width = 6.8, height = 5, units = "in")
perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

filtered \(\rho = 0.50\)

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.50) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

filtered \(\rho = 0.90\)

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.90) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

filtered \(\rho = 0.05\)

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = model)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size),
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = rel_bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RPB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

perf_crit  |> filter(var_combo_graph == 0.05) |> 
  ggplot(aes(x = variance, y = bias_icc, fill = nj_size)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(n_bar_prop ~ n_bar_size, 
             labeller = label_bquote(rows = zeta ==.(n_bar_prop), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "PB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

3 SE of ICC

3.1 Interactions

3.1.1 True ICC vs num pooled ICCs

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)),
                                     cols = k ==.(icc_est_n))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       # strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RSEB_numICC.png", width = 6.8, height = 5, units = "in")
# perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
#   geom_hline(yintercept = .05, linetype = "dashed")  + 
#   geom_hline(yintercept = -.05, linetype = "dashed")  + 
#   geom_boxplot(alpha = .5) + 
#   facet_grid( var_combo_graph_f ~ icc_est_n, 
#              labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
#                                             cols = k ==.(icc_est_n))
#              ) + 
#   labs(x = NULL, y = "RSEB", fill = NULL) + 
#   theme_bw() +
#   theme(
#     legend.position = "bottom",
#         axis.text.x = element_text(size = 6.5),
#         axis.text.y = element_text(size = 6.5),
# 
#         plot.caption=element_text(hjust = 0, size = 6),
#        # legend.title=element_text(size = 6), 
#         legend.text=element_text(size = 6.5),
#         axis.title.x = element_text(size=6.5),
#         axis.title.y = element_text(size=6.5),
#       #  strip.background =element_rect(fill="white")
#         strip.text.x = element_text(size = 6.5, face="bold"),
#         strip.text.y = element_text(size = 6.5, face="bold")
#       )

#ggsave("RPB_se_numICC.png", width = 6.8, height = 5, units = "in")
perf_crit |>
  ggplot(aes(x = as.character(variance), y = rel_bias_se_ICC, 
             color = model, 
             shape = as.character(tau)
             )
         ) +
  geom_point(alpha = .8, 
           #  position = position_jitter(height = .3), 
             size = 1)   + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  +
  facet_grid( var_combo_graph_f ~ icc_est_n, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = k ==.(icc_est_n))) + 
  labs(
    x = NULL, 
    y ="RSEB",
    color = "", shape = ""
  ) + 
  theme_bw() +
  theme(
    plot.caption=element_text(hjust = 0, size = 10),
    legend.position= "bottom",
    panel.spacing.x = unit(5, "mm"), 
    panel.spacing.y = unit(5, "mm"),
  )

ggsave("RSEB_numICC_alt2.png", width = 6.8, height = 5, units = "in")

3.1.2 True ICC vs tau (between study )

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ tau, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = tau==.(tau))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RSEB_tau.png", width = 6.8, height = 5, units = "in")

3.1.3 True ICC vs number of clusters

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ nj_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = j== .(nj_size))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RPB_se_j.png", width = 6.8, height = 5, units = "in")

3.1.4 True ICC vs cluster size

TASHA – interestingly for SE of the pooled estimate, despite the bias in the pooled estimate.

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_size, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RSEB_nj.png", width = 6.8, height = 5, units = "in")

3.1.5 True ICC vs Degree of unbalance

perf_crit  |> 
  ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( var_combo_graph_f ~ n_bar_prop, 
             labeller = label_bquote(rows = rho ==.(as.character(var_combo_graph_f)), 
                                            cols = zeta ==.(n_bar_prop))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RSEB_n_bar_prop.png", width = 6.8, height = 5, units = "in")

3.1.6 number clusters vs cluster size

perf_crit  |> ggplot(aes(x = variance, y = rel_bias_se_ICC, fill = model)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid( nj_size ~ n_bar_size, 
             labeller = label_bquote(rows = j ==.(nj_size), 
                                            cols = bar(n_j)==.(n_bar_size))) + 
  labs(x = NULL, y = "RSEB", fill = NULL) + 
  theme_bw() +
  theme(
    legend.position = "bottom",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),

        plot.caption=element_text(hjust = 0, size = 6),
       # legend.title=element_text(size = 6), 
        legend.text=element_text(size = 6.5),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
      #  strip.background =element_rect(fill="white")
        strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

#ggsave("RSEB_j_nj.png", width = 6.8, height = 5, units = "in")

4 L2 Variance

4.1 Main Effects

4.1.1 L2 variance

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
 # facet_grid( rows = vars(var_combo_graph_f)
              #labeller = label_parsed
#              ) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6, face="bold"),
        plot.caption=element_text(hjust = 0, size = 8),
        legend.title=element_text(size = 8), 
        legend.text=element_text(size = 8),
        axis.title.x = element_text(size=8),
        axis.title.y = element_text(size=8),
        strip.background =element_rect(fill="white"))

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus_ratio)) + 
  geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
 # facet_grid( rows = vars(var_combo_graph_f)
              #labeller = label_parsed
#              ) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6, face="bold"),
        plot.caption=element_text(hjust = 0, size = 8),
        legend.title=element_text(size = 8), 
        legend.text=element_text(size = 8),
        axis.title.x = element_text(size=8),
        axis.title.y = element_text(size=8),
        strip.background =element_rect(fill="white"))

4.1.2 L2 variance by tau

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(tau )
             , 
             labeller = label_bquote(cols = tau ==.(tau))) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_L2_tau_mean.png", width = 6.8, height = 5, units = "in")



perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus_ratio)) + 
  geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(tau )
             , 
             labeller = label_bquote(cols = tau ==.(tau))) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

4.1.3 L2 variance by number of clusters

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
#  geom_hline(yintercept = 0, linetype = "solid", colour = "red")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(nj_size )
             , 
             labeller = label_bquote(cols = j ==.(nj_size))) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_L2_j_mean.png", width = 6.8, height = 5, units = "in")



perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_clus_ratio)) + 
  geom_hline(yintercept = 1, linetype = "solid")  + 
 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(nj_size )
             , 
             labeller = label_bquote(cols = j ==.(nj_size))) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

4.1.4 L1 variance by clustersize

perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_res)) + 
  geom_hline(yintercept = .05, linetype = "dashed")  + 
  geom_hline(yintercept = -.05, linetype = "dashed")  + 
#  geom_hline(yintercept = 0, linetype = "solid", colour = "red")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(n_bar_size )
             , 
             labeller = label_bquote(cols = nj ==.(n_bar_size))) + 
  labs(x = "ICC Parameter", y = "RPB") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

ggsave("RPB_L2_nj_mean.png", width = 6.8, height = 5, units = "in")



perf_crit2  |> ggplot(aes(x = factor(var_combo_graph), y = rel_bias_res_ratio)) + 
   geom_hline(yintercept = 1, linetype = "solid")  + 
  geom_boxplot(alpha = .5) + 
  facet_grid(cols = vars(n_bar_size )
             , 
             labeller = label_bquote(cols = nj ==.(n_bar_size))) + 
  labs(x = "ICC Parameter", y = "RPB Ratio") + 
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 6.5),
        axis.text.y = element_text(size = 6.5),
      #  plot.caption=element_text(hjust = 0, size = 8),
      #  legend.title=element_text(size = 8), 
      #  legend.text=element_text(size = 8),
        axis.title.x = element_text(size=6.5),
        axis.title.y = element_text(size=6.5),
       strip.text.x = element_text(size = 6.5, face="bold"),
        strip.text.y = element_text(size = 6.5, face="bold")
      )

5 Evaluating the Distribution ICC estimates

Note for the Fisher TF pooled estimates, they have already been back transformed for all of these analyses.

load("data/results_conv.RData")
load("data/results_pooled_est.RData")

6 Distributions of the ICC Estimates from primary studies

skewness <- function(x) {
  
  variance = var(x)
  K = length(x)
  skewness = (1/(K*sqrt(variance)^3)*sum((x-mean(x))^3))
  
  return(skewness)
  
}

kurtosis <- function(x) {
  
  variance = var(x)
  K = length(x)
  kurtosis = (1/(K*variance^2))*sum((x-mean(x))^4)
  
  return(kurtosis)
  
}

First, which set of conditions resulted in the lowest raw bias and relative bias?

For \(\rho = 0.5\), it is when \(k = 50\), \(n_j \sim U[30, 50]\), \(\overline{n_j} \sim U[10, 30]\), \(n_{prop} = 0.1\), and \(\tau = 0.2\). I will use that set for the distributions.

# perf_crit |>  group_by(ICC_true) |>  slice(which.min(bias_icc)) |> select(ICC_true, method, var_icc_name, icc_est_n, nj_size, n_bar_size, n_bar_prop, tau)
# 
# perf_crit |>  group_by(ICC_true) |>  slice(which.min(rel_bias_icc)) |> select(ICC_true, method, var_icc_name, icc_est_n, nj_size, n_bar_size, n_bar_prop, tau)


results_all_converged_dist_icc_est <- results_all_converged |> filter(icc_est_n == 50) |>  filter(nj_size == "U[30, 50]") |>  filter(n_bar_size == "U[10, 30]") |>  filter(n_bar_prop == 0.1) |>  filter(tau == 0.02) 


# results_all_converged_dist_icc_est <- results_all_converged |> filter(nj_size == "U[30, 50]") |> filter(n_bar_size == "U[10, 30]") |> 
#   filter(n_bar_prop == 0.1) |> filter(tau == 0.02 ) |> filter(icc_est_n == 50) 

6.1 Distribution of ICC estimates when \(\rho = 0.05\):

dist_0.05 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.05) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.05  <- data.frame(icc = dist_0.05)

ggplot(dist_0.05 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.05, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.05$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.05$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .05 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.0500555, variance of 0.0010352, median of 0.0464974, a skewness of 0.5933149, and kurtosis of 3.1094325.

6.2 Distribution of ICC estimates when \(\rho = 0.10\):

dist_0.10 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.10) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.10  <- data.frame(icc = dist_0.10)

ggplot(dist_0.10 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
      geom_vline(xintercept = 0.1, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.10$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.10$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .10 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.0993899, variance of 0.0015262, median of 0.0964493, a skewness of 0.4208593, and kurtosis of 3.1433641.

6.3 Distribution of ICC estimates when \(\rho = 0.15\):

dist_0.15 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.15) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.15  <- data.frame(icc = dist_0.15)

ggplot(dist_0.15 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
      geom_vline(xintercept = 0.15, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.15$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.15$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .15 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.1482727, variance of 0.0020559, median of 0.1457706, a skewness of 0.3036363, and kurtosis of 3.0933228.

6.4 Distribution of ICC estimates when \(\rho = 0.25\):

dist_0.25 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.25) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.25  <- data.frame(icc = dist_0.25 )

ggplot(dist_0.25 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.25, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.25$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.25$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .25 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.2471755, variance of 0.0031014, median of 0.245689, a skewness of 0.140802, and kurtosis of 2.9323949.

6.5 Distribution of ICC estimates when \(\rho = 0.50\):

dist_0.5 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.5) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.5  <- data.frame(icc = dist_0.5 )

ggplot(dist_0.5 , aes(x = icc)) +
 geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.50, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.5$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.5$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .50 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.4923617, variance of 0.0044361, median of 0.4947604, a skewness of -0.2104669, and kurtosis of 3.0122678.

compare this to a distribution that has the larger average cluster sizes condition:

results_all_converged_dist_icc_est2 <- results_all_converged |> filter(nj_size == "U[30, 50]") |> filter(n_bar_size == "U[30, 50]") |> 
  filter(n_bar_prop == 0.1) |> filter(tau == 0.02   ) |> filter(icc_est_n == 50) 

dist_0.5 <- results_all_converged_dist_icc_est2 |> ungroup() |>  filter(var_combo_graph == 0.50) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.5  <- data.frame(icc = dist_0.5)

ggplot(dist_0.5 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.50, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.5$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.5$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .50 and cluster size 'U[30, 50]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.4928036, variance of 0.0040805, median of 0.4950266, a skewness of -0.1914525, and kurtosis of 2.9920759.

6.6 Distribution of ICC estimates when \(\rho = 0.90\):

dist_0.9 <- results_all_converged_dist_icc_est |> ungroup() |>  filter(var_combo_graph == 0.9) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.9  <- data.frame(icc = dist_0.9)

ggplot(dist_0.9 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.9, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.9$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.9$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .9 and cluster size 'U[10, 30]'", x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

This distribution has a mean of 0.8956405, variance of 0.0011112, median of 0.8993834, a skewness of -0.691453, and kurtosis of 3.6654658.

dist_0.9 <- results_all_converged_dist_icc_est2 |> ungroup() |>  filter(var_combo_graph == 0.9) |> select(clustervar, residvar) |> distinct()  |>
  mutate(icc_list = map2(clustervar, residvar, ~ map2(.x, .y, ~ .x / (.x + .y)))) |>  select(icc_list) |> unnest(icc_list) |>  unlist() 

dist_0.9  <- data.frame(icc = dist_0.9)

ggplot(dist_0.9 , aes(x = icc)) +
  geom_histogram(fill = "#92D895", color = "#e9ecef", alpha = 0.6, binwidth = 0.01) +
  geom_vline(xintercept = 0.9, linetype = "solid") + 
  geom_vline(xintercept = mean(dist_0.9$icc), linetype = "dashed") + 
  geom_vline(xintercept = median(dist_0.9$icc), linetype = "dotted") + 
  labs(title = "Histogram of ICC Values for rho = .9 and cluster size 'U[30, 50]'" , x = "ICC Estimate", y = "Frequency") +
  theme_minimal()

7 Distributions of the Pooled ICC estimates

7.1 \(\rho = 0.05\):

results_all_converged |> filter(var_combo_graph == 0.05) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.05,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.05)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")
variance model K Mean Pooled Estimate var_est Raw Bias Pooled Estimate MCSE Raw Bias RMSE min_est max_est
Donner REML 48000 0.0444794 1.93e-05 -0.0055206 2.00e-05 0.0070531 0.0194512 0.0689376
Donner RVE 48000 0.0444223 1.95e-05 -0.0055777 2.02e-05 0.0071169 0.0197097 0.0696369
Fisher REML 48000 0.0444866 1.93e-05 -0.0055134 2.00e-05 0.0070468 0.0197216 0.0690194
Fisher RVE 48000 0.0444274 1.95e-05 -0.0055726 2.02e-05 0.0071105 0.0198642 0.0696210
Fisher TF REML 48000 0.0461986 1.67e-05 -0.0038014 1.87e-05 0.0055821 0.0241912 0.0702339
Fisher TF RVE 48000 0.0461986 1.67e-05 -0.0038014 1.87e-05 0.0055822 0.0241917 0.0702424
Hedges REML 48000 0.0445867 1.92e-05 -0.0054133 2.00e-05 0.0069609 0.0196622 0.0691946
Hedges RVE 48000 0.0445332 1.94e-05 -0.0054668 2.01e-05 0.0070212 0.0199373 0.0698602
Smith REML 48000 0.0443622 1.94e-05 -0.0056378 2.01e-05 0.0071521 0.0195930 0.0687574
Smith RVE 48000 0.0443009 1.96e-05 -0.0056991 2.02e-05 0.0072169 0.0196059 0.0693541
Swiger REML 48000 0.0444893 1.93e-05 -0.0055107 2.00e-05 0.0070442 0.0197166 0.0690243
Swiger RVE 48000 0.0444302 1.95e-05 -0.0055698 2.02e-05 0.0071078 0.0198623 0.0696227
results_all_converged |> filter(var_combo_graph == 0.05) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance ) 

7.2 \(\rho = 0.10\):

results_all_converged |> filter(var_combo_graph == 0.10) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.10,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.10)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est), 
            .groups = "drop")|> 
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")
variance model K Mean Pooled Estimate var_est Raw Bias Pooled Estimate MCSE Raw Bias RMSE min_est max_est
Donner REML 48000 0.0929020 2.93e-05 -0.0070980 2.47e-05 0.0089287 0.0638537 0.1205014
Donner RVE 48000 0.0929014 2.93e-05 -0.0070986 2.47e-05 0.0089284 0.0637351 0.1203858
Fisher REML 48000 0.0928594 2.95e-05 -0.0071406 2.48e-05 0.0089687 0.0636333 0.1204919
Fisher RVE 48000 0.0928520 2.95e-05 -0.0071480 2.48e-05 0.0089749 0.0635221 0.1203793
Fisher TF REML 48000 0.0950363 2.59e-05 -0.0049637 2.32e-05 0.0071069 0.0673182 0.1239894
Fisher TF RVE 48000 0.0950365 2.59e-05 -0.0049635 2.32e-05 0.0071069 0.0673183 0.1239837
Hedges REML 48000 0.0930215 2.93e-05 -0.0069785 2.47e-05 0.0088291 0.0640675 0.1205711
Hedges RVE 48000 0.0930298 2.92e-05 -0.0069702 2.47e-05 0.0088208 0.0639626 0.1204625
Smith REML 48000 0.0926788 2.97e-05 -0.0073212 2.49e-05 0.0091247 0.0636260 0.1202392
Smith RVE 48000 0.0926590 2.97e-05 -0.0073410 2.49e-05 0.0091417 0.0635153 0.1200836
Swiger REML 48000 0.0928633 2.94e-05 -0.0071367 2.48e-05 0.0089651 0.0636440 0.1204934
Swiger RVE 48000 0.0928562 2.94e-05 -0.0071438 2.48e-05 0.0089710 0.0635328 0.1203808
results_all_converged |> filter(var_combo_graph == 0.1) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance )

7.3 \(\rho = 0.15\):

results_all_converged |> filter(var_combo_graph == 0.15) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.15,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.15)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|> 
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")
variance model K Mean Pooled Estimate var_est Raw Bias Pooled Estimate MCSE Raw Bias RMSE min_est max_est
Donner REML 48000 0.1418869 4.15e-05 -0.0081131 2.94e-05 0.0103601 0.1064030 0.1766836
Donner RVE 48000 0.1418123 4.17e-05 -0.0081877 2.95e-05 0.0104256 0.1063915 0.1766836
Fisher REML 48000 0.1418120 4.17e-05 -0.0081880 2.95e-05 0.0104287 0.1063188 0.1766414
Fisher RVE 48000 0.1417274 4.19e-05 -0.0082726 2.95e-05 0.0105041 0.1063003 0.1766414
Fisher TF REML 48000 0.1436102 3.84e-05 -0.0063898 2.83e-05 0.0089018 0.1106306 0.1773362
Fisher TF RVE 48000 0.1436106 3.84e-05 -0.0063894 2.83e-05 0.0089014 0.1106319 0.1773445
Hedges REML 48000 0.1420110 4.14e-05 -0.0079890 2.94e-05 0.0102592 0.1068142 0.1768486
Hedges RVE 48000 0.1419499 4.15e-05 -0.0080501 2.94e-05 0.0103125 0.1068147 0.1768486
Smith REML 48000 0.1416008 4.20e-05 -0.0083992 2.96e-05 0.0106104 0.1056585 0.1765680
Smith RVE 48000 0.1414930 4.23e-05 -0.0085070 2.97e-05 0.0107085 0.1056643 0.1765680
Swiger REML 48000 0.1418164 4.17e-05 -0.0081836 2.95e-05 0.0104245 0.1063331 0.1766487
Swiger RVE 48000 0.1417324 4.19e-05 -0.0082676 2.95e-05 0.0104994 0.1063149 0.1766487
results_all_converged |> filter(var_combo_graph == 0.15) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model~ variance )

7.4 \(\rho = 0.25\):

results_all_converged |> filter(var_combo_graph == 0.25) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.25,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.25)^2) / K),
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>

  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")
variance model K Mean Pooled Estimate var_est Raw Bias Pooled Estimate MCSE Raw Bias RMSE min_est max_est
Donner REML 48000 0.2413431 6.46e-05 -0.0086569 3.67e-05 0.0118127 0.2013939 0.2871551
Donner RVE 48000 0.2411995 6.50e-05 -0.0088005 3.68e-05 0.0119352 0.2016610 0.2870117
Fisher REML 48000 0.2412389 6.49e-05 -0.0087611 3.68e-05 0.0119028 0.2012907 0.2869823
Fisher RVE 48000 0.2410835 6.54e-05 -0.0089165 3.69e-05 0.0120371 0.2013082 0.2868213
Fisher TF REML 48000 0.2410603 6.57e-05 -0.0089397 3.70e-05 0.0120662 0.1997265 0.2822733
Fisher TF RVE 48000 0.2410612 6.57e-05 -0.0089388 3.70e-05 0.0120653 0.1997266 0.2822763
Hedges REML 48000 0.2414566 6.45e-05 -0.0085434 3.67e-05 0.0117264 0.2011547 0.2875039
Hedges RVE 48000 0.2413259 6.49e-05 -0.0086741 3.68e-05 0.0118374 0.2017192 0.2873750
Smith REML 48000 0.2410222 6.54e-05 -0.0089778 3.69e-05 0.0120820 0.2012777 0.2869650
Smith RVE 48000 0.2408420 6.60e-05 -0.0091580 3.71e-05 0.0122404 0.2012916 0.2868021
Swiger REML 48000 0.2412436 6.49e-05 -0.0087564 3.68e-05 0.0118986 0.2012878 0.2869940
Swiger RVE 48000 0.2410888 6.54e-05 -0.0089112 3.69e-05 0.0120324 0.2013215 0.2868341
results_all_converged |> filter(var_combo_graph == 0.25) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model~ variance )

7.5 \(\rho = 0.50\):

results_all_converged |> filter(var_combo_graph == 0.50) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
      #      "Median Pooled Estimate" = median(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.50,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.50)^2) / K), 
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")
variance model K Mean Pooled Estimate var_est Raw Bias Pooled Estimate MCSE Raw Bias RMSE min_est max_est
Donner REML 48000 0.4953113 0.0000912 -0.0046887 4.36e-05 0.0106376 0.4411449 0.5419561
Donner RVE 48000 0.4953182 0.0000913 -0.0046818 4.36e-05 0.0106394 0.4411045 0.5420179
Fisher REML 48000 0.4952857 0.0000913 -0.0047143 4.36e-05 0.0106546 0.4410808 0.5419028
Fisher RVE 48000 0.4952923 0.0000914 -0.0047077 4.36e-05 0.0106570 0.4410339 0.5419671
Fisher TF REML 48000 0.4880025 0.0001072 -0.0119975 4.73e-05 0.0158473 0.4341028 0.5367767
Fisher TF RVE 48000 0.4880040 0.0001072 -0.0119960 4.73e-05 0.0158458 0.4340954 0.5367603
Hedges REML 48000 0.4953049 0.0000911 -0.0046951 4.36e-05 0.0106353 0.4412048 0.5418731
Hedges RVE 48000 0.4953117 0.0000912 -0.0046883 4.36e-05 0.0106370 0.4411663 0.5419321
Smith REML 48000 0.4952682 0.0000915 -0.0047318 4.37e-05 0.0106725 0.4411312 0.5421544
Smith RVE 48000 0.4952744 0.0000916 -0.0047256 4.37e-05 0.0106756 0.4409836 0.5422489
Swiger REML 48000 0.4952861 0.0000913 -0.0047139 4.36e-05 0.0106541 0.4410833 0.5419024
Swiger RVE 48000 0.4952927 0.0000914 -0.0047073 4.36e-05 0.0106565 0.4410366 0.5419665
results_all_converged |> filter(var_combo_graph == 0.5) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance )

7.6 \(\rho = 0.90\):

results_all_converged |> filter(var_combo_graph == 0.90) |> 
  group_by(variance, model) |> 
  summarise(K = n(),
            "Mean Pooled Estimate" = mean(pooled_icc_est),
      #      "Median Pooled Estimate" = median(pooled_icc_est),
            var_est = var(pooled_icc_est),
            "Raw Bias Pooled Estimate" = mean(pooled_icc_est) - 0.90,
            "MCSE Raw Bias" = sqrt(var_est / K),
            RMSE = sqrt(sum((pooled_icc_est - 0.90)^2) / K), 
            min_est = min(pooled_icc_est), 
            max_est = max(pooled_icc_est),
            .groups = "drop")|>
  kbl() |>
  kable_styling() |> 
  scroll_box(width = "900px", height = "500px")
variance model K Mean Pooled Estimate var_est Raw Bias Pooled Estimate MCSE Raw Bias RMSE min_est max_est
Donner REML 48000 0.9026873 1.98e-05 0.0026873 2.03e-05 0.0051952 0.8735445 0.9278472
Donner RVE 48000 0.9026360 1.98e-05 0.0026360 2.03e-05 0.0051716 0.8732500 0.9277722
Fisher REML 48000 0.9027883 1.99e-05 0.0027883 2.03e-05 0.0052565 0.8735703 0.9280164
Fisher RVE 48000 0.9027390 1.99e-05 0.0027390 2.03e-05 0.0052311 0.8732787 0.9279004
Fisher TF REML 48000 0.8970148 2.40e-05 -0.0029852 2.24e-05 0.0057375 0.8673503 0.9201010
Fisher TF RVE 48000 0.8970153 2.40e-05 -0.0029847 2.24e-05 0.0057371 0.8673525 0.9201015
Hedges REML 48000 0.9025850 1.98e-05 0.0025850 2.03e-05 0.0051422 0.8735030 0.9274451
Hedges RVE 48000 0.9025347 1.98e-05 0.0025347 2.03e-05 0.0051211 0.8732056 0.9274583
Smith REML 48000 0.9029853 2.00e-05 0.0029853 2.04e-05 0.0053753 0.8738152 0.9280318
Smith RVE 48000 0.9029397 2.00e-05 0.0029397 2.04e-05 0.0053475 0.8735355 0.9279116
Swiger REML 48000 0.9027842 1.99e-05 0.0027842 2.03e-05 0.0052539 0.8735697 0.9280039
Swiger RVE 48000 0.9027347 1.99e-05 0.0027347 2.03e-05 0.0052285 0.8732780 0.9278907
results_all_converged |> filter(var_combo_graph == 0.9) |> 
#  group_by(var_combo_graph) |> 
  ggplot(aes(x = pooled_icc_est)) +
  geom_histogram( fill = "#404080", color = "#e9ecef", alpha = 0.6, position = 'identity', bins = 50) +
    facet_grid(model ~ variance )

7.7 Skewness and Kurtosis

perf_crit |> 
  group_by(variance, model, var_combo_graph) |> 
  summarise(m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  pivot_wider(names_from = variance, values_from = c(m_s, m_k)) |>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value",
                    "Donner",
                    "Fisher",
                    "Fisher TF",
                    "Hedges",
                    "Smith",
                    "Swiger",
                    "Donner",
                    "Fisher",
                    "Fisher TF",
                    "Hedges",
                    "Smith",
                    "Swiger"
                    )) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "Mean Skewness" = 6, "Mean Kurtosis" = 6)) |> 
  column_spec(9:14, bold = T, background = "grey") 
Mean Skewness
Mean Kurtosis
Model ICC Value Donner Fisher Fisher TF Hedges Smith Swiger Donner Fisher Fisher TF Hedges Smith Swiger
REML 0.05 0.0592022 0.0594254 0.0866294 0.0593649 0.0591826 0.0594198 2.954853 2.955471 2.953104 2.954562 2.955909 2.955497
REML 0.10 0.0723233 0.0724675 0.0749330 0.0718898 0.0727667 0.0724434 3.001821 3.001935 2.995341 3.001267 3.002176 3.001900
REML 0.15 0.0357151 0.0359730 0.0377721 0.0360568 0.0364047 0.0359696 2.947726 2.948172 2.955624 2.947198 2.948294 2.948137
REML 0.25 0.0323212 0.0323009 0.0280801 0.0319869 0.0327775 0.0322975 3.010076 3.009571 3.014416 3.008855 3.010905 3.009555
REML 0.50 -0.0362609 -0.0361881 -0.0318783 -0.0357684 -0.0362934 -0.0361812 2.980498 2.980601 2.978749 2.980647 2.980997 2.980596
REML 0.90 -0.0973374 -0.0970725 -0.1078941 -0.0969041 -0.0974283 -0.0970904 3.027450 3.027156 3.022494 3.027109 3.026623 3.027172
RVE 0.05 0.0684156 0.0684702 0.0866546 0.0682543 0.0689004 0.0684497 2.951688 2.952241 2.953011 2.951607 2.952394 2.952227
RVE 0.10 0.0742096 0.0746966 0.0749193 0.0736681 0.0752857 0.0746794 2.998217 2.998166 2.995317 2.997296 2.998940 2.998164
RVE 0.15 0.0359919 0.0364623 0.0378290 0.0359824 0.0374692 0.0364526 2.949006 2.949565 2.955679 2.948395 2.950002 2.949530
RVE 0.25 0.0336099 0.0338578 0.0281223 0.0329711 0.0346274 0.0338303 3.005622 3.004860 3.014341 3.005028 3.005913 3.004849
RVE 0.50 -0.0364144 -0.0363523 -0.0318745 -0.0358920 -0.0365080 -0.0363443 2.980415 2.980558 2.978801 2.980482 2.981113 2.980552
RVE 0.90 -0.1068451 -0.1071497 -0.1078788 -0.1061320 -0.1078974 -0.1071320 3.024866 3.024611 3.022530 3.024555 3.023916 3.024628
# %>%
#   row_spec(3:5, bold = T, color = "white", background = "#D7261E")


perf_crit |> 
  filter(var_combo_graph == 0.05) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.05" = 4))
ICC = 0.05
Model ICC Value Mean Bias Mean Relative Bias Mean Skewness Mean Kurtosis
REML Donner -0.0055206 -0.1104120 0.0592022 2.954853
REML Fisher -0.0055134 -0.1102687 0.0594254 2.955471
REML Fisher TF
-0.003801
TF
-0.003801
perf_crit |> 
  filter(var_combo_graph == 0.1) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.10" = 4))
ICC = 0.10
Model ICC Value Mean Bias Mean Relative Bias Mean Skewness Mean Kurtosis
REML Donner -0.0070980 -0.0709803 0.0723233 3.001821
REML Fisher -0.0071406 -0.0714056 0.0724675 3.001935
REML Fisher TF
-0.004963
TF
-0.004963
perf_crit |> 
  filter(var_combo_graph == 0.15) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.15" = 4))
ICC = 0.15
Model ICC Value Mean Bias Mean Relative Bias Mean Skewness Mean Kurtosis
REML Donner -0.0081131 -0.0540873 0.0357151 2.947726
REML Fisher -0.0081880 -0.0545870 0.0359730 2.948172
REML Fisher TF
-0.006389
TF
-0.006389
perf_crit |> 
  filter(var_combo_graph == 0.25) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.25" = 4))
ICC = 0.25
Model ICC Value Mean Bias Mean Relative Bias Mean Skewness Mean Kurtosis
REML Donner -0.0086569 -0.0346275 0.0323212 3.010076
REML Fisher -0.0087611 -0.0350445 0.0323009 3.009571
REML Fisher TF
-0.008939
TF
-0.008938
perf_crit |> 
  filter(var_combo_graph == 0.5) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.50" = 4))
ICC = 0.50
Model ICC Value Mean Bias Mean Relative Bias Mean Skewness Mean Kurtosis
REML Donner -0.0046887 -0.0093774 -0.0362609 2.980498
REML Fisher -0.0047143 -0.0094286 -0.0361881 2.980601
REML Fisher TF
-0.011997
TF
-0.011996
perf_crit |> 
  filter(var_combo_graph == 0.9) |>
  group_by(model, variance) |> 
  summarise(
            m_b = mean(bias_icc),
            m_rb = mean(rel_bias_icc),
            m_s = mean(skewness),
          #  min_s = min(skewness),
          #  max_s = max(skewness),
            m_k = mean(kurtosis),
          #  min_k = min(kurtosis),
          #  max_k = max(kurtosis)
          .groups = "drop")|>
  kable(booktabs = TRUE, col.names = c("Model",
                    "ICC Value", "Mean Bias", "Mean Relative Bias", "Mean Skewness", "Mean Kurtosis")) |>
  kable_styling() |> 
  add_header_above(c(" " = 2, "ICC = 0.90" = 4))
ICC = 0.90
Model ICC Value Mean Bias Mean Relative Bias Mean Skewness Mean Kurtosis
REML Donner 0.0026873 0.0029859 -0.0973374 3.027450
REML Fisher 0.0027883 0.0030982 -0.0970725 3.027156
REML Fisher TF
-0.002985
TF
-0.002984
perf_crit_agg <- perf_crit |> 
  group_by(model, variance, icc_est_n , var_combo_graph, tau) |> 
  summarise(skewness_mean = mean(skewness), 
            kurtosis_mean = mean(kurtosis),
            n = n(),
            .groups = "drop")



ID1 <- c(0.05, 0.10, 0.15)

perf_crit_agg |> filter(var_combo_graph < 0.2) |> ggplot(aes( var_combo_graph, skewness_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 0 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID1), breaks = ID1) +
  theme_minimal()

ID2 <- c(0.25, 0.5, 0.9)

perf_crit_agg |> filter(var_combo_graph > 0.2) |> ggplot(aes( var_combo_graph, skewness_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 0 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID2), breaks = ID2) +
  theme_minimal()

ID1 <- c(0.05, 0.10, 0.15)

perf_crit_agg |> filter(var_combo_graph < 0.2) |> ggplot(aes( var_combo_graph, kurtosis_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 3 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID1), breaks = ID1) +
  theme_minimal()

ID2 <- c(0.25, 0.5, 0.9)

perf_crit_agg |> filter(var_combo_graph > 0.2) |> ggplot(aes( var_combo_graph, kurtosis_mean, col=variance, shape = model ) ) + 
  facet_grid(tau ~ icc_est_n, 
             labeller = label_bquote(rows = tau == .(tau), 
                                            cols = k ==.(icc_est_n))) +

  geom_point( alpha=0.75 ) + 
  geom_line( alpha=0.75 ) +
  geom_hline( yintercept = 3 ) + 
  scale_x_continuous("ICC Parameter", labels = as.character(ID2), breaks = ID2) +
  theme_minimal()

Other Graphs

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = factor(icc_est_n))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) + 
  geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ model, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = .(model))) +
  labs(x = "ICC Value", y = "Skewness", color = "k") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = factor(tau))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Skewness", color = "tau") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = nj_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Skewness", color = "") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = skewness, color = n_bar_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
  geom_hline(yintercept = 0, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Skewness", color = "") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = factor(icc_est_n))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
  geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ model) + 
  labs(x = "ICC Value", y = "Kurtosis", color = "k") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = factor(tau))) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Kurtosis", color = "tau") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = nj_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
    geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Kurtosis", color = "") 

perf_crit |> 
  ggplot(aes(x = factor(var_combo_graph), y = kurtosis, color = n_bar_graph)) +
  geom_point(alpha = .5, position = position_jitter(width = .3)) +
  geom_hline(yintercept = 3, linetype = "solid")  + 
  facet_grid(variance ~ icc_est_n, 
             labeller = label_bquote(rows = .(variance), 
                                            cols = k ==.(icc_est_n))) +
  labs(x = "ICC Value", y = "Kurtosis", color = "")